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Abstract 

We present ah initio calculations of the high-temperature axial c/a ratio of hexagonal- 
close-packed (hep) iron at Earth's core pressures, in order to help interpret the ob- 
served seismic anisotropy of the inner core. The calculations are based on density 
functional theory, which is known to predict the properties of high-pressure iron with 
good accuracy. The temperature dependence of c/a is determined by minimising the 
Helmholtz free energy at fixed volume and temperature, with thermal contributions 
due to lattice vibrations calculated using harmonic theory. Anharmonic corrections 
to the harmonic predictions are estimated from calculations of the thermal average 
stress obtained from ah initio molecular dynamics simulations of hep iron at the 
conditions of the inner core. We find a very gradual increase of axial ratio with 
temperature. This increase is much smaller than found in earlier calculations, but 
is in reasonable agreement with recent high-pressure, high-temperature diffraction 
measurements. This result casts doubt on an earlier interpretation of the seismic 
anisotropy of the inner core. 
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1 Introduction 



In the last few years, there have been a number of ah initio studies of iron 
at high temperatures and pressures, in an effort to understand the proper- 
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ties of the Earth's solid inner c ore (IStixrude et a 



(119951) :ISoderlind et al. 
IVocadl o et all (Il997l): 



Wasserman et alJ 



1(ll994l) : IStixrude and CohenI 
(ll996h:IStixrude et all (ligm 
1999); IVocadlo et al. 



■■_ Steinle-Neumann et al.l ^..^^ 

Belonoshko et all (120031) :lLaio et all (|2nnnl ): lAlfe et all f.2001'1: Steinle-Neum ann et al. 



(|200lL |2002[ ): IVocadlo et al.1 (|2003[ )). including its elastic properties. A contro- 
versial issue, addr essed by some o f this w o rk, co ncerns the elastic anisotropy 
of the inner core (Creager (|l992h : lTroml (|l993t) ). i.e. the fact that compres- 
sional waves are observed to traverse the core region some 3-4% faster along 
the Earth's rotational axis than in the equatorial plane. It is widely assumed 
that the phase of Fe in the inner core is hexagonal close packed (hep), and it 
has become clear that to understand the elastic properties it is necessary to 
know how the axial ratio c/a depends on temperature and pressure. In order 
to provide information about this, we present here ab initio calculations of the 
axial ratio of hep iron over a range of pressures and temperatures relevant to 
the inner core. 



Ah initio calculations on pure crystalline iron are very relevant for under- 
standing the inner core, even though the core is known to contain light im- 
purities (the lead ing can d idates are O, Si and S), a nd is believed to c ontain 
Ni as well as Fe (|Poirieil (|l994l) ). Recent estimates (|Alfe et all (|2002al) ) sug- 
gest that the fraction of light impurities in the inner core is approximately 
8.5 mol%, and this low fraction will probably not change its properties greatly. 
The Ni content will also probably make only a small difference, since the 
electronic structures of Ni and Fe are so similar. The determination of the 
properties of pure Fe is thus essential as a starting point for further refine- 
ments. The common assumption tha t Fe in the inn e r core i s in the hep struc- 



ture h as recently been challenged (jVocadlo et al. (|2003t ): iBelonoshko et al 
( 2003[ )). and it is possible that the crys t al str u cture is different in d ifferent 



parts of the core (ISong and Helmbergei ( 1998[ ): Ishii and Dziewonski (Eb02 ). 
Beghein and Tramperd ( 2003[ )). However, in order to make progress, it is clearly 
essential to understand the properties of the leading candidates, of which hep 
i s certainly one. It is well es t ablished that density f u nctional theory (^DFT) 
(iHohen berg and Kohn ( 1964 ): Kohn and Sham ( 1965[ l: Jones and Gunnarsson 
(|1989) ) gives a rather accurate description of Fe and other transition metals, 
both under ambient conditions, and under high pressures and temperatures. 
This is known from comparisons with experiment for a wide range of prop- 
erties, including the equilibrium lattice parameter under ambient conditions, 
elastic constants, mag n etic properties, phonon frequencies (IStixrude et al. ( 1994 ) ; 
Soderhnd et~aD (|l99fil) ; IVocadlo et all (|l997l) ; IVocadlo et al.l (|2000>) ): thepressure- 
volume re lationship at high pressu r es, the phonon density of sta t es an d the 
Hugoniot (|Wasserman et al.l (|l99fi[ ): iMao et al.1 (|200l[ ): lAlfe et all (|200lh . 



A number of groups have used DFT calculations to investigate the elastic con- 
stants of hep Fe under inne r-core conditions. The early c alculations were per- 
formed at zero temperature ( Stixrude and CohenI |l995.) : .Steinle-Neumann et al. 
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and led to the suggestion that elastic anisotropy could be explained by 
a preferential alignment of crystalline c-axes with the Earth's rotational axis. 
However, more recent work has indicated a st rong temperature dependence of 
the elastic constants (ISteinle-Neumann et al.l f(20QL.2002 )). including a cross- 
ing of the cii and C33 constants at approximately 1500 K, which would result 
in a reversal of the necessary crystalline alignment. It appears that the strong 
temperature dependence of the elastic constants is due in part to a rapid 
increase of the axial ratio with increasing temperature. However, this work 
can be questioned, because it is based upon a stati stical mechanical a pprox - 
i mation known as th e particle- in-cell (PIC) model (iHirshfelder et al. (Eosi ): 



iHolt and Ros?j (ll970l] : |Holt et al.l()l97i , ^, 

f)l97.5h : ICowlev et all (|l990l )). which may not be reliable. In a ddition, indepen 



Ree and Holtl (|l973h : IWestra and Cowlev 



dent P IC calculations have failed to reprod uce their results (jCannarelli et al 
Very recently, experimental work ()Ma et al.l (|2004[ l) has shown that 



this strong c/a variation with temperature is not observed at pressures up to 
160 GPa. Our aim in this paper is to perform calculations of the high temper- 
ature axial ratio, which cirG , clS far as possible, free of statistical mechanical 
approximat ions . 



The calculation of phonon dispersion relations using DFT is nowadays com- 
pletely routine, and this makes it possible to obtain the Helmholtz free energy 
F of a high-temperature crystal without any statistical mechanical approxi- 
mations, provided that anharmonicity can be ignored. Phonon dispersion re- 
lations for Fe in both the magnetic bo dy-centred cub i c and the hep structures 
have been reported by several groups (IVocadlo et all feOOOh ). The equilibrium 
value of c/a can then be obtained for given V and T by minimising F with 
respect to c/a. Of course, at temperatures near the melting point, it is not 
obvious that anharmonic corrections to F can be neglected. However, direct 
DFT molecular dyn amics simulations al low the calculation of the thermal av- 
erage stress tensor ( Oganov et al. ( 200l[ )). so that the equilibrium c/a can be 
found by requiring that the stress tensor be isotropic; such direct simulations 
completely include anharmonicity. The DFT calculation of phonon frequencies 
is straightforward and relatively inexpensive, but ah initio molecular dynam- 
ics demands much larger computational resources. Our strategy in this work 
is therefore to base most of the calculations on the harmonic approximation, 
but to use a small number of molecular dynamics simulations to estimate 
anharmonic corrections to the equilibrium c/a. 



The rest of this paper is organised as follows: In Sec. 2 we present the technical 
background to our electronic structure calculations, and describe the statistical 
mechanical principles of our harmonic and molecular-dynamics calculations. 
In Sec. 3 we describe our harmonic calculations, including the technical tests 
we performed to ensure their accuracy. In Sec. 4 we present our ah initio 
molecular-dynamics results. Discussion of our results for the temperature and 
pressure dependence of the axial ratio follows in Sec. 5. 
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2 Techniques 



The principal quantity in this work is the ab initio Helmholtz free energy 
F{V, q, T) as a function of the volume V , axial ratio q = c/a and temperature 
T . Since the temperatures of interest here are far above the Debye tempera- 
ture, F is given by the standard expression of classical statistical mechanics 



F(V,g,T) = -fceT'ln j dri . . . drN Qw[-Uxi{ji . . .Yn]T^\) /k^T]^ (1) 

where T is the temperature and A is the thermal wavelength. The quantity 
?7Ai(ri . . . r^v; Tel) is the ah initio Helmholtz free energy of the system when the 
atoms are fixed at the positions {rj}; U^i is a free energy because it includes 
thermal excitation of the electrons, at electronic te mperatu r e Tpi, as described 
by finite temperature density functional theory ( Mermin (W6^). Since we 
are interested in the case of full thermal equilibrium, we have T^i = T. In 
the present work, the DFT calculation of Uai is made using the exchange- 
correlation energy given by the PW91 generalised gradient approximation 
fjWang and Perdewl (|l99lh l. The imp l ementation of DFT is the p r oiecto r aug- 
mented wave (PAW) scheme fiBlochll ( 1994 ): Kresse and Joubert ('1999')') with 
core radii, augmentation charge radii, etc. as reported in (|Alfe et all (j20^). 
As in our previous work, atomic states up to and including 3p are treated as 
core states, but the high-pressure response of 3s and 3p states is included via 
an empirical pair potential. This approximation is tested using spot-checks 
with 3p electrons inclu ded in the valance set. A l l calcu lations were performed 
using the VASP code ( Kresse and Furthmiiller ()l996bl lah). 



In practice, the ab initio Helmholtz free energy is separated into perfect lattice, 
harmonic and anharmonic contributions: 



anhami 



(2) 



Here -Fperf = ?7ai(Ri . . . Rat; Tgi) is the ab initio free energy with all atoms 
fixed on the ir perfect-lattice p ositions {Ri}. The harmonic contribution -Fharm 



is given by ()Alfe et all (|200lO) 



Fharm — S/CbTIu 



with u the geometric mean phonon frequency: 

1 



In a; 



(3) 



(4) 



where iVks is the total number of /c-points and phonon-branches, the sum being 
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taken over the first Brillouin zone. is a function of V, q and Te\. As usual, the 
phonon frequencies uj]^s at each wavevector k are obtained by diagonahsing the 
dynamical matrix, whose elements are defined in terms of the force-constant 
matrix ^isa,i't(3 as: 



1 

M 



lsa,l'tf3 exp 



(5) 



where M is the mass of each atom. Here is the equilibrium position of 
the sth atom in the Zth primitive cell. Calculation of th e force-consta nt ma- 
trix $ is perforra e d by t he srn all- displacement method (|Kresse et al ] (|W95); 
Alfe et all (|2nni[) : IXlfj (|l998[ )). in which each atom in the primitive cell is 
given a small displacement, and DFT is used to calculate the resulting force 
on every atom in a large repeating cell. The calculation of the dynamical 
matrix and hence the mean frequency Co must be taken to convergence with 
respect to the size of the repeating cell. For the hep structure, the entire 
$ matrix is cal culated by perfor ming only two independent displacements. 



as explained in lAlfe et al 



()2001f ). Within the harmonic approximation, the 
equilibrium c/a at given atomic volume and temperature is determined by 
calculating Fperf(g) and -Fharm(9) at a series of q values, and by minimising the 
quantity Fperf(g) + Fharm(g) numerically. 

Our estimates of anharmonic corrections to the equilibrium c/a are based on ab 
initio molecular dynamics simulations. T hese are pe r forme d in the canonical 
ensemble, using an Andersen thermostat (Andersen ( 198Cl| )). The determina- 
tion of the equilibrium c/ a is based on a calculation of the time-averaged stress 
tensor cTq,/?. The system is in hydrostatic equilibrium with respect to variation 
of q when (733 — (Xn = 0. In comparing with harmonic predictions, it is useful 
to note that the components of the thermal average stress tensor are related 
to the ab initio harmonic free energy by 



— lim 



dF 



de 



0/3 / np 



(6) 



where Sap is the infinitesimal strain tensor and Vq is the volume of the sys- 
tem before application of the strain. In particular, for our constant volume 
calculations of F{y, q, T), we can write 



< 0-33 - 0-11 >= — — 



3q (OF' 
W \ 'dq^ 



(7) 



VT 



This means that the stress component < (T33 — an > within the harmonic 
approximation can be obtained by taking the derivative with respect to q of 
the free energy Fperf(g) + Fharm(g)- 
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3 Harmonic calculation of the equilibrium axial ratio 



In order to determine the equilibrium axial ratio in harmonic theory, the two 
quantities that it is necessary to calculate are the perfect lattice free energy 
Fperf, and the geometric mean phonon frequency both as function of volume, 
axial ratio c/a and electronic temperature. The calculation of these quantities 
is presented in the following subsections, and the results for the equilibrium 
axial ratio are presented in Sec. 3.3. The precision we need to achieve in 
the calculations is set by the precision with which we wish to determine the 
equilibrium axial ratio. As a guideline, we set ourselves the target of obtaining 
c/a to a precision of ±0.005. 



3. 1 Perfect lattice free energy 



DFT r esults for the per fect lattice free energy for c/a = 1.6 were reported ear- 



l ier bylAlfe et al.l (|2001l l. and were more recently reported bv lGannarelli et al. 
( 20031 ) for q in the range 1.48-1.72. The present calculations are closely re- 
lated to the previous ones, but we have introduced some refinements. We have 
calculated -Fperf for electronic temperatures in the range 1000-7000 K. For T^x 
in this range, we find that an 8 x 8 x 5 Monkhorst-Pack fc-point sampling set 
provides a precision of better than 1.5 meV atom^^. The plane- wave cutoff has 
been set to achieve complete continuity in the curve Fperf(g), and to ensure 
that at every point on this curve, the calculated energy and stress components 
are converged to within our required targets. 



For the set of Tei mentioned above, we have calculated Fperf for c/a in the 
range 1.54-1.70, at steps of 0.01 and for volumes between 6.8 and 8.8 at 
steps of 0.2 A^, as well as at three particular volumes of interest, 6.97, 7.50 and 
8.67 A^. A parameterisation essentially exactly fitting the data consists of a 
cubic polynomial in c/a, whose coefficients are fitted to quartic polynomials in 
Tel. As a point of reference for our later discussion, we have used these results 
to calculate the equilibrium axial ratio at T = as a function of volume. We 
have converted the results to obta i n q as a function of pressure by using the 
P{y) relation from lStixrude et al.l (|l994^ . In Fig. 1, we com pare the r e sultin g 
values of q{P) with very recent synchrotron measurements ( Ma et al. ( 20041) ) 
in the pressure range 60-160 GPa. In both cases, we see a very slight increase 
of q with pressure, though the theoretical values are lower by about 0.008. 
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Fig. 1. Calculated equilibrium axial ratio of hep Fe at zero tei nperature as a fu nction 
of pressure (solid curve) , compared with diffraction r esults of lMa et al.l ^2004 ) (open 
squares with error bars). 



3.2 Harmonic calculations 



As explained in Sec. 2, the harmonic vibrational free energy is determined 
entirely by u. Values of u must be converged with respect to A;-point sampling, 
cell size and atomic displacement, and we must ensure that this convergence 
is achieved across the required range of c/a. In order to achieve the required 
tolerance in the harmonic free energy, systematic errors in ui must not vary 
with respect to q by more than 0.1%. 

We have designed a strategy for ensuring convergence of u. In our initial 
tests, we ado pt a displacemen t of 0.01 A, which has been used in previous 
calculations 



) pt a gispiacemen t oi u.Ui A, wnicn nas Deen usea m previous 
(|Alfe et all (|200lh . With this displacement, we determine the 



A;-point sampling needed to obtain convergence for a 16-atom supercell. We 
then study convergence with respect to cell size, the fc-point sampling for 
each cell size being chosen in the light of the /c-point sampling tests of the 
16-atom system. This set of tests tells us the cell size needed. To conclude 
the tests, we then address convergence of atomic displacement and fc-point 
sampling for this cell size. All these tests are performed at T^i = 1000 K. This 
temperature is deliberately chosen to be much lower than the temperatures of 
real interest. We do this because the fineness of A;-point sampling needed to 
achieve a given degree of convergence decreases with increasing temperature. 
One can therefore be sure that if convergence is achieved with respect to k- 
point sampling at 1000 K, it will certainly be achieved at higher temperatures. 
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For the 16-atom system, Gj was calculated for Monkhorst Pack sets up to 
9x9x6. We find that for a Monkhorst Pack 5x5x3 grid, Uj is converged to 
within 0.01%. Next, a set of supercells was chosen, such that the dimension 
along the the hexagonal axis is similar to the dimension in the basal plane. 
Sizes chosen were 16 (2 x 2 x 2), 36 (3 x 3 x 2), 54 (3 x 3 x 3), 96 (4 x 4 x 3), 128 
(4x4x4) and 150 (5x5x4) atoms. For each cell size, we calculate Uj using 
fc-point sampling that is equivalently finer than the converged value of 5 x 5 x 3 
in the 16-atom system. We find that for a 54-atom system, uj is converged to 
within 0.5%, however the consistency between calculations for different values 
of c/a at this cell size, is much better than our tolerance. Finally, fc-point and 
displacement convergence were carried out on the 54-atom system. We found 
that for a 3 X 3 X 2 Monkhorst Pack grid and a displacement of 0.01 A, uj is 
converged to within approximately 1% in absolute terms, but, again, to well 
within our tolerances for non-cancelling errors. These tests were all performed 
for both c/a = 1.60 and 1.70. All the harmonic results that follow were based 
upon values of oj calculated for these values of supercell size, fc-point sampling 
and displacement. 

Calculations of uj were performed for atomic volumes of 6.97, 7.50 and 8.97 A^, 
for c/a from 1.58 to 1.70 in steps of approximately 0.03, and for electronic 
temperatures of 2000, 4000 and 6000 K. We parameterised lnci;(g) for each 
volume and temperature using a second order polynomial in g, which gave an 
effectively exact fit to the results. 



3. 3 Harmonic results for equilibrium axial ratio 



Equilibrium values of c/a were obtained at 6.97, 7.50 and 8.67 A^ and tem- 
peratures of 0, 2000, 4000 and 6000 K, by analytic minimisation of the total 
harmonic free energy Fperf(g) + -Fharm(Q'), using the polynomial parameterisa- 
tions described above for In Lj(y, q, T^i) and Fperf (V, q, T^i). Results are shown in 
Fig. 2 , together with the previous theoretical resul ts of ISteinle- Neumann et al 



( 2001 ) and the very recent experimental results of iMa et al. ( 20041 ). Note that 
all the theoretical results show the variation of c/a with T at fixed volume, 
whereas the experimental results are at the fixed pressure of 161 GPa. The 
present results differ greatly from the earlier theoretical results, in that we 
find only a very moderate increase of c/a with temperature. The experimental 
increase of c/a with te r apera ture is also far smaller than the predictions of 



Steinle-Neumann et al.l (|2001l l. and is, if anything, smaller than the variation 



that we predict. More detailed discussion will be given in Sec. 5. 
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Fig. 2. Calculated equilibrium axial ratio as a function of temperature for different 
volumes. For this work (light curves) atomic volumes are 6.97 (solid curve 



7.50 (dashed curve) and 8.67 A^ (dotted curve). For Steinle-Neumann et al 



( 200lh (heavy curves) volumes are 6.81 A^ (solid curve), 7.11 A^ (dashed c urve) and 
7.41 A ^ (dotted curve). Also shown are diffraction measurements due to iMa et al 



(|2Qi]J) at 7.73 A^/atom (open squares with error bars) and our own ab initio MD 
calculations at 6.97 A^. 

4 Molecular dynamics calculations 



The determination of the equilibrium axial ratio by ab initio molecular dynam- 
ics is based on the calculation of stresses; specifically on the stress difference 
o"33 — cii, which disappears for the equilibrium value of c/a. In order to apply 
this technique successfully, a number of technical sources of error must be 
brought under control. Firstly, the duration of the run must be long enough 
to reduce statistical errors to within an acceptable level; secondly, the size 
of the simulation cell must be large enough; and thirdly, electronic /c-point 
sampling must be adequate. The first two of the convergence questions can 
be addressed in detail using classical molecular dynamics simulations, and 
we have performed tests of this kind as described below. The question of k- 
point convergence will be discussed when we present our ab initio molecular 
dynamics calculations in Sec. 4.2. 

Although ab initio molecular dynamics calculations fully include anharmonic- 
ity, they are very demanding in terms of computational resources, so that we 
can only perform these calculations at a very small number of state points. 
Our aim is to put limits on the possible size of anharmonic corrections to the 
stress, and hence to c/a. Our ab initio molecular dynamics results for a^^ — au 
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can be compared with our harmonic calculations by using Eqn. (7) to calculate 
C33 — c"ii in the harmonic approximation. Since anharmonicity is negligible at 
low temperatures, the harmonic and molecular dynamics calculations of the 
stress should be in close agreement under these conditions. This provides a 
further check both on our harmonic and molecular dynamics results. 



4.1 Molecular dynamics simulations for an embedded- atom model 



Our tests on statistical and system-size errors were perfor med using an em- 



bedde d atom model (EAM) for Fe, similar to that used bv iBelonoshko et al 



n the EAM, the total energy ^'tot of a system of N atoms is rep- 
resented as the sum of two contributions. First, an empirical, repulsive pair 
potential, and secondly, a sum of embedding energies for each atom, repre- 
senting quantum mechanical band-structure effects: 



1 



pair 2 



-E'embed = f{Pi) 



Pi=Yl ^(ki-ril) 
where 0(r) and ipir) have the inverse power forms 



r) = e - ) and 



(9) 



and f{pi) = —eC{pi)2. In practice, both and ip are cut off smoothly by 
making the transformation 



r) + (3 + •yr r < Tq 
3 



a[ri — r 




rQ < r < ri , 
ri < r 



(10) 



where ro and ri are cutoff radii, with a, (3 and 7 chosen such that 0, cf)' and 
(f)" are continuous at tq. The form guarantees that these three functions are 
also continuous at ri. In the present calculations, we use ri = 7 A, tq = 
0.9ri. The pararneters have been modified somewhat from those given by 
Belonoshko et al. ( 2000l ). in order to reproduc e better ou r ah initio mo lecular 
dynamics simulations on solid and liquid Fe (|Alfe et al.1 (|200lL l2002bl )). The 
parameters we user are: n = 5.93, m = 4.788, e = 0.1662 eV, a = 3.4714 A 
and C = 16.55. 
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30 100 150 

Number of atoms 

Fig. 3. 15 ps timc-avcraged calculations of ao = (733 — an for a variety of cell sizes 
within the embedded atom model. 

In order to assess the length of run needed to reduce the statistical error on 
— o"!! to an acceptable level, we have performed molecular dynamics runs at 
several system sizes and at a number of different thermodynamic state points. 
An illustration of such a run is given in Fig. 4. On the basis of these calculations 
we find that, in order to reduce the r.m.s. statistical errors to within 1 GPa, for 
example, we require a run of length 3-4 ps. Since the r.m.s. error is inversely 
proportional to the square root of the length of run, we may immediately infer 
the length of run necessary to obtain any required precision. From statistical 
mechanics, we expect that for a given length of run, the r.m.s. error on the 
stress components will decrease with increasing system size as N~2^ where N 
is the number of atoms. We find that this is consistent with tests performed 
on system size, with r.m.s. errors on a 1.5 ps average falling from 2 GPa for a 
54-atom cell, to 1 GPa for 150 atoms. To test size errors we performed runs on 
systems containing 54, 96, 128 and 150 atoms at identical state points. Fig. 3 
shows the effect of system size on (J33 — an. Error bars show r.m.s. statistical 
errors as described above. We see that for a system of 96 atoms, the error due 
to size effects is within approximately 0.6 GPa. 

It follows from these tests that in order to calculate 0-33 — an to a tolerance 
needed to determine the equilibrium axial ratio to within 0.005, it should be 
enough to perform ab initio molecular dynamics simulations of 1 ps on 96 
atoms. 

As a further test of our techniques, we have compared molecular dynamics 
results for (T33 — an with the predictions of harmonic theory for the embedded 
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10 20 30 

Time (ps) 



Fig. 4. Instantaneous (dotted curve) and average (solid) value of ctq = (T33 — an 
(upper plot) and standard error (lower plot) for a 40 ps run for V = 6.97 AVatom, 
T = 4000 K and c/a = 1.65 in the embedded atom model. 

atom model. Fig. 5 presents the molecular dynamics and harmonic results for 
V = 6.97 A^, c/a = 1.65 for a range of temperatures. Because of the gradual 
variation of ass — an with temperature, we use longer runs of up to 5 ps. Cell 
sizes are converged withing the limits discussed above. We see that at low 
temperatures, the molecular dynamics results reproduce well the predictions 
of harmonic theory. This agreement provides additional confirmation that the 
statistical and system-size errors on the molecular dynamics calculations are 
very small. The deviation, rising to around 8.5 GPa at 6000 K, shows the 
high-temperature emergence of anharmonic effects. 



4-2 Ah initio molecular dynamics results 

We have performed ah initio molecular dynamics simulations at an atomic 
volume of 6.97 A^, in which wc have calculated (733 — an. The molecular 
dynamics calculations arc performed using exactly the same DFT methods as 
in the harmonic calculations. From the tests on the embedded atom model, 
we already have a good indication of the cell size and length of run necessary, 
but the question of /c-points must also be addressed. To do this, we calculate 
o'ss — <^ii for several disordered configurations, selected from classical M.D. 
trajectories, both with F-point sampling and with larger numbers of /c-points. 
Our tests are performed for a 96-atom system, with an atomic volume of 
6.97 A^ and an axial ratio of 1.65. The molecular dynamics simulation from 
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Fig. 5. Molecular dynamics (open circles) and harmonic (diamonds) calculations of 
= 1^33 ~ (^11 for 6.97 cell in the embedded atom model, with c/a = 1.70. 
The dashed curve represents a second-order fit to the MD data, while the solid line 
is its tangent at OK. In agreement with thermodynamic theory, this coincides the 
harmonic calculations of the stress. 



which the configurations were drawn was performed at 5000 K. We took 4 
snapshots at intervals of 0.2 ps. For this state point, (733 — (Th, with F-point 

sampling is equal to 7.0 ± 2 GPa. For a 2 x 2 x 2 /c-point grid, it is equal to 
7.1 ± 2.6 GPa. Similar tests were performed to ensure that size effects were 
consistent with the embedded atom model. 

Our ab initio calculations at 6.97 were performed for c/a — 1.615 and 
1.66, and for T = 4000 K and 5000 K. Fig. 6 shows results for (733 — an 
along with harmonic predictions. We see that our molecular dynamics results 
closely resemble the predictions of harmonic theory, and vary in the same way 
with c/a. However, even allowing for the error bars, there is an appreciable 
difference between the harmonic and M.D. results. One of the technical sources 
of error which could account for this discrepancy is system size effects. In order 
to check this, we have repeated the 5000 K runs using 150 atoms instead of 96. 
The results of this test suggest that system-size errors could account for most 
of this discrepancy. We have also examined other technical sources of error, 
such as plane-wave cutoff and /c-point sampling; however, as described above, 
there are not capable of producing such a discrepancy. It is possible in principle 
that some of this discrepancy can be attributed to anharmonic effects, but it 
is difficult to separate this from other effects with any certainty. Even allowing 
for the remaining uncertainties, it seems certain that anharmonic effects are 
not capable of shifting the equilibrium axial ratio by more than about 0.01. 
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1.62 



c/a 



1.64 



1.66 



Fig. 6. Comparison of ab initio M.D. (heavy curves) and harmonic (light) results for 
the stress (T33 — fxn. Solid lines are at 4000 K, dashed lines at 5000 K. The dotted line 
shows results for 150-atom molecular dynamics at 4000 K. The diamond represents 
a 96-atom run including 3p states explicitly in the valence set. All calculations 
were performed for an atomic volume of 6.97 A^ . For comparison with earlier PIC 
calculations due to ISteinle-Neumann et al.l (|200ll ). see Fig. 2. 



We have also used ab initio M.D. to test another significant question. All our 
harmonic and M.D. calculations up to this point have treated 3p and 3s elec- 
trons as core states, but with the high-pressure response of these states treated 
via an empirical pair potential. In order to test the effect of this approximation, 
we have repeated the 96-atom M.D. calculations at 5000 K with 3p electrons 
explicitly included in the valence set. The resulting stress — an agrees with 
our other results within our error bars. This indicates that the effects of this 
approximation are negligible for our present purposes. 



5 Conclusions 



We have shown that it is rather straightforward to calculate harmonic free 
energies, and hence calculate the variation of c/a within the harmonic ap- 
proximation. Because these calculations are computationally inexpensive, it is 
possible to apply them to a wide variety of thermodynamic state points. Since 
our earlier comparisons with experiment ( Mao et al. ( 200l[ )) show the reliabil- 
ity of our density functional methods for calculating the phonon frequencies 
for iron at high pressures, we expect our predictions of the temperature depen- 
dence of c/a to be reliable. We have also shown that direct ab initio molecular 
dynamics simulation can be used to put limits on the possible size of anhar- 
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monic effects. These methods are rather expensive, but provide a means to 
perform "spot checks" on our results. The molecular dynamics calculations 
show that anharmonic corrections to the c/a ratio are small, even at temper- 
atures near the melting point. 

The main scientific conclusion from the work is that c/a varies only rather 
weakly with temperature, and that its variation becomes weaker as pressures 
approach those of the inner core. Quantitatively, at a pressure of 100 GPa, c/a 
increases from 1.585 at zero temperature, to 1.61 at the melting point, and at 
300 GPa, from 1.59 at zero temperature, to 1. 62 at the melting point. These 



results do not support the earlier predictions bv lSteinle-Neumann et al.l (|200l| ) 
of a strong increase of c/a with temperature within the particle-in-cell approx- 
imation. Interestingly, the results we have found ar e very similar to our ow n 



predictions from the particle-in-cell approximation ( Gannarelli et al. ( 2003[ )) 



This implies that this approximation cannot be respo nsible for the large dis- 
crepan cy between the present results and those of Steinle-Neumann et al.l 



( 20011 1. This finding of a rather weak increase of c/a wit h temperature i s 



strongly supported by the recent diffraction experiments of iMa et al.l (|2004 ). 
Even allowing for the significant error bars on their c/a results, it seems clear 
that at 161 GPa, 2000 K, c/a is no bigger than about 1.61 at most. Our results 
would indicate a value of around 1.60 under these conditions. 

Our result s have implications for understanding the elastic anisotropy of the 
inner core. Steinle-Neumann et al.l ((2001.) predicted a reversal of the crystalline 



alignment required to explain the anisotropy of the inner core, due to a crossing 
of the Cii and C33 elastic moduli at approximately 1500 K. They attribute this 
to the strong temperature-dependence they find in c/a. Since our results cast 
doubt on this strong temperature dependence, they also cast doubt on the 
proposed explanation for the elastic anisotropy. However, in order to resolve 
this question fully, we need to make ah initio predictions of the elastic moduli, 
that are free of statistical mechanical approximations. The methods we have 
presented here should be capable of achieving this, and we hope to be able to 
present such results in due course. 
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